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Analysis: Modeling Event Times by a 
Stochastic Process Reaching a Boundary 
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Abstract. Many researchers have investigated first hitting times as 
models for survival data. First hitting times arise naturally in many 
types of stochastic processes, ranging from Wiener processes to Markov 
chains. In a survival context, the state of the underlying process repre- 
sents the strength of an item or the health of an individual. The item 
fails or the individual experiences a clinical endpoint when the process 
reaches an adverse threshold state for the first time. The time scale can 
be calendar time or some other operational measure of degradation or 
disease progression. In many applications, the process is latent (i.e., 
unobservable) . Threshold regression refers to first-hitting-time models 
with regression structures that accommodate covariate data. The pa- 
rameters of the process, threshold state and time scale may depend on 
the covariates. This paper reviews aspects of this topic and discusses 
fruitful avenues for future research. 
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ing risk, cure rate, duration, environmental studies, first hitting time, 
gamma process, lifetime, latent variable models, maximum likelihood, 
operational time, occupational exposure, Ornstein-Uhlenbeck process, 
Poisson process, running time, stochastic process, stopping time, sur- 
vival analysis, threshold regression, time-to-event, Wiener diffusion pro- 
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1. INTRODUCTION 

Many types of lifetime, duration or time-to-event 
data may be interpreted as first hitting times (FHT's) 
of a boundary or threshold state by sample paths of 
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a stochastic process, which may be latent or observ- 
able. FHT models have a long history of application 
in diverse fields, including medicine, environmental 
science, engineering, business, economics and soci- 
ology. FHT models may describe the length of a 
hospital stay, the survival time of a transplant pa- 
tient, the onset time for a cancer induced by occu- 
pational exposure, the failure time of an engineering 
system, the depletion time of an inventory, the sur- 
vival time of a new business, the transition time for 
a stock price change and the length of a marriage. 
Relevant articles, covering theory and application, 
appear in both the lifetime data and reliability lit- 
eratures. These FHT models are gradually becom- 
ing widely adopted because of their conceptual ap- 
peal, realism and applicability. Recently, interest in 
them has deepened and spread, and exciting new 
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areas of application are being encountered. The po- 
tential applications require new conceptual view- 
points, theoretical advances, analytical techniques 
and methodological extensions to which the discus- 
sion returns later. 

To make FHT models truly valuable in applica- 
tions, they must be capable of extension to include 
regression structures. Regression structures allow the 
effects of covariates to explain the inherent disper- 
sion of the data, thereby taking account of vari- 
ability and sharpening inferences. Regression struc- 
tures also provide scientific insights into potential 
causal roles of covariates in the underlying processes, 
boundary sets and time scales. This article reviews 
aspects of FHT models and is concerned especially 
with regression structures for FHT models, which 
will be referred to as threshold regression, or TR 
for short. The word "threshold" refers to the fact 
that the FHT is triggered by the underlying process 
reaching a threshold state within a boundary set, as 
described in more detail in the next section. 

2. THE FIRST-HITTING-TIME (FHT) MODEL 

A first-hitting-time (FHT) model has two basic 
components: (1) a parent stochastic process {X(t),t G 
T,x G X} with initial value X(0) = xq, where T is 
the time space and X is the state space of the pro- 
cess and (2) a boundary set B, where B C X . We 
shall refer to the boundary set B boundary, 
barrier or threshold, depending on which is most de- 
scriptive or conventional in the context. The process 
{X(t)} may have a variety of properties, such as one 
or many dimensions, the Markov property, continu- 
ous or discrete states, and monotonic sample paths. 
Whether the sample path of the parent process is 
observable or latent (unobservable) is an important 
distinguishing characteristic of the FHT model. La- 
tent processes are the most common by far. The 
boundary set B may also have different features as 
will be illustrated in later examples. 

Taking the initial value X(0) = xq of the process 
to lie outside the boundary set B, the first hitting 
time of B is the random variable S defined as 

(1) S = m£{t:X(t) eB}. 

Thus, the first hitting time is the time when the 
stochastic process first encounters set B. We refer 
to the state first encountered in the boundary set 
by the process, that is, X(S) S B, as the threshold 
state. The boundary set defines a stopping condition 



for the process and, therefore, the FHT is usually 
a stopping time in the formal sense of that term in 
stochastic process theory. Note that when the parent 
process is latent, there is no direct way of observing 
the FHT event in the state space of the process. 

In some versions of the FHT model, there is no 
guarantee that process {X(t)} will reach the bound- 
ary set B, so P(S < oo) < 1. We will let S = oo de- 
note the absence of a finite hitting time with P(S = 
oo) = 1 — P(S < oo). The later discussion will show 
situations where this condition is plausible and a 
desirable model feature. The basic FHT model (1) 
assumes that B is fixed in time. In some applica- 
tions, however, it varies with time, that is, B(t). This 
variation may be deterministic or follow a stochastic 
process. 

An exhaustive review of the first-hitting-time lit- 
erature is impossible within the confines of a single 
article. Eaton and Whitmore (1977) discuss FHT's 
general model for hospital stay. Aalen and 
Gjessing (2001) provide an excellent overview of much 
of this subject. Likewise, Lawless (2003) gives a com- 
plete and compact summary of theory, models and 
methods (see Section 11.5, pages 518-523). Lee and 
Whitmore (2004) also provide an overview of first- 
hitting-time models for survival and time-to-event 
data. We will make numerous references to selected 
work as we proceed. There is a huge literature deal- 
ing with theoretical and mathematical aspects of 
FHT models that we will not attempt to review or 
incorporate. We also will not cover random growth 
curve models, such as those of Carey and Koenig 
(1991) and Lu and Meeker (1993), which have an 
FHT interpretation but where the only randomness 
at the level of the individual parent process is con- 
fined to a noise factor. The large literature on linear 
and nonlinear regression methods for survival data 
and reliability where the underlying models have no 
central FHT motivation (such as accelerated failure 
time and proportional hazards models) is also not 
covered. 

3. EXAMPLES OF FIRST-HITTING-TIME 
MODELS 

The parent stochastic processes may take many 
forms, from Wiener processes to Markov chains. Like- 
wise, the nature of the boundary state may vary 
widely — for example, a fixed threshold in a Wiener 
process or an absorbing state in a Markov chain. 
As the preceding description of a first-hitting-time 



THRESHOLD REGRESSION MODELS 



3 



model is quite abstract, we now list a few basic ex- 
amples to illustrate the variety encountered in ap- 
plications. 

1. Bernoulli process and negative binomial first 
hitting time. The number of trials S required to 
reach the mth success in a Bernoulli process {Bt,t = 
1,2,...} has a negative binomial distribution with 
parameters m and p, where p is the success prob- 
ability on each trial. To give this setup our stan- 
dard representation, we consider a parent process 
{Xt, t = 0, 1, 2, . . .} with initial value Xq = xq = m 
and let X t = x$ — B t ,t = 1,2, ... , where {B t } is the 
preceding Bernoulli process. The first hitting time 
is the first Bernoulli trial t = S for which Xt equals 
0. The number of rocket launches required to get m 
satellites in orbit is a simple example of this FHT 
model. 

2. Poisson process and Erlang first hitting time. 
The time S until the occurrence of the mth event in 
a Poisson process {N(t),t > 0} with rate parameter 
A has an Erlang distribution with parameters m and 
A. Again, to give this setup our standard represen- 
tation, we consider a parent process {X(t),t > 0} 
with initial value X(0) = xq = m and let X(t) = 
xo — N(t), where {N(t),t > 0} is the preceding Pois- 
son process. The first hitting time is the earliest time 
t = S when X(t) = 0. This FHT model is illustrated 
by the time to failure of an engineering system con- 
sisting of m components in parallel, having identi- 
cal and independent exponential lifetimes, that are 
placed in service successively as failures occur. 

3. Wiener process and inverse Gaussian first hit- 
ting time. Consider a Wiener process {X(t),t > 0} 
with mean parameter /i, variance parameter a 2 , and 
initial value X(0) = xq > 0. The time S required for 
the process to reach the zero level for the first time 
has an inverse Gaussian distribution if the process 
mean parameter [i is negative so the process tends 
to drift toward the zero level. Lancaster (1972) de- 
scribes the duration of an industrial strike using 
this model. In this context, {X(t)} represents the 
distance between the positions of management and 
labor at time t after the start of the strike. The 
initial separation of the parties is X(0) = xq > 0. 
The strike ends when the process first encounters 
the zero level where the parties agree and settle. In 
addition to Lancaster (1972), refer to applications 
described in Whitmore and Neufeldt (1970), Whit- 
more (1975, 1979, 1983, 1986, 1995), Doksum and 
Hoyland (1992), Doksum and Normand (1995), Lu 



(1995), Whitmore and Schenkelberg (1997) and Hor- 
rocks and Thompson (2004), to name a few. Onar 
and Padgett (2000) and Padgett and Tomlinson (2004) 
extend the Wiener diffusion model to an acceler- 
ated testing context. Pettit and Young (1999) set 
the model in a Bayesian context. As illustrated later, 
the inverse Gaussian distribution has a closed-form 
probability density function (p.d.f.) and a compu- 
tationally simple cumulative distribution function 
(c.d.f.). Their formulas vary slightly depending on 
whether the parent process is defined as rising or 
falling to hit the relevant boundary. 

4. Gamma process and inverse gamma first hitting 
time. Consider a parent process {X{t),t > 0} with 
initial value X(0) = xq > 0. Let X(t) = xq — Z{t) 
where {Z(t),t > 0} is a gamma process with scale 
parameter (3, shape parameter a and Z(0) = 0. The 
first hitting time of the zero level in the parent pro- 
cess (X = 0) has an inverse gamma c.d.f., defined 
by the identity P(S > t) = P{Z{t) < x ). The iden- 
tity follows from the fact that a gamma process 
has monotonic (nondecreasing) sample paths. The 
availability of computational routines for the gamma 
c.d.f. allows the c.d.f. of S to be computed read- 
ily. Singpurwalla (1995) and Lawless and Crowder 
(2004) consider the gamma process as a model for 
degradation. Park and Padgett (2005) consider both 
geometric Brownian motion and gamma processes in 
an accelerated degradation model. 

5. Ornstein-Uhlenbeck process and Ricciardi-Sato 
first hitting time. The Ornstein-Uhlenbeck (OU) pro- 
cess is a variant of a Wiener process that is mean- 
reverting in that it tends to drift back toward a fixed 
equilibrium level and thus has a property of home- 
ostasis. Aalen and Gjessing (2004) propose the first- 
hitting-time distribution for such a process as a sur- 
vival model and derive pertinent results. They point 
out that the form of the FHT distribution is found 
in Ricciardi and Sato (1988), who have studied it 
extensively. 

6. Markov chain and absorbing state first hitting 
time. Markov chains {Xt, t = 0, 1, 2, . . .} are an im- 
portant type of parent process. The state space X 
consists of the possible states of the chain. The time 
space T is the transition steps for the chain. The 
first hitting time is the minimum number of steps 
required to move from an initial state Xo = xo to 
a set of boundary states B. The FHT distribution 
depends on the transition matrix of the chain in a 
precise mathematical manner. As a case example, a 
Markov chain can model product brand switching in 
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the field of consumer behavior. An FHT of interest 
might be the number of purchases that will be made 
in the product category before a consumer who cur- 
rently uses brand xq will switch to another brand, 
say, b. In this case, B is the singleton set {b}. 

7. Semi-Markov processes and their first hitting 
times. All of the preceding examples are special cases 
of Markov processes. A semi-Markov process {X(t), 
t > 0} extends the Markov chain model by includ- 
ing the random time that the process resides in each 
state. Although the Markov property is generally 
lost by this extension, the model remains of great 
practical value. In a semi-Markov model, the first 
hitting time represents the time that the process 
resides in the initial and subsequent states before 
it first enters one of the states that define set B. 
There are many important examples of this multi- 
state model, some more complicated than others. 
The two-stage clonal conversion model for cancer 
provides a case example. This is a chemical carcino- 
genesis model that, in its basic form, postulates a 
pool of (stem) cells that are susceptible to malig- 
nant transformation. The cells proceed through an 
initiation stage and then through malignant conver- 
sion according to a two-stage Markov process with 
fixed transition rates. Once a malignant cell is gen- 
erated, a second statistical model describes progres- 
sion through stages of active cancer to death (home- 
ostasis, angiogenesis, metastasis, death). The model 
has been elaborated in a series of works by Mool- 
gavkar, Luebeck and Anderson (1998) and Luebeck 
et al. (1999), and others. The model has been ap- 
plied, for example, in a study of lung cancer risk 
posed to Chinese tin miners by arsenic, radon and 
tobacco exposure (Hazelton et al., 2001). 

4. LATENT FIRST HITTING TIMES AND 
COMPETING RISKS 

Most duration data are gathered under conditions 
of competing risks in which two or more causes are 
competing to determine the observed duration. The 
outcome becomes associated with both a time and 
mode of occurrence. For example, a medical first- 
hitting-time model that describes the time of a sub- 
ject's death may recognize explicitly that multiple 
causes are competing to produce death. Both the 
time and cause of death are recorded for each sub- 
ject. Refer to Kalbfleisch and Prentice (1980, 2002) 
and Crowder (2001), for example, for detailed tech- 
nical discussions of this topic. 



FHT models accommodate the competing risk as- 
pect in a natural fashion. In an FHT model with 
competing risks, the boundary set B is partitioned 
into mutually exclusive and exhaustive subsets, say 
Si, ... , Be, associated with FHT causes c = 1, . . . , C. 
Let D denote the observed cause. Then realization 
D = d associated with the observed FHT outcome 
S would be defined as 

(2) D = d if X(S)eB d . 

The concept of a latent FHT offers an interesting 
vehicle for discussing competing risks. In this frame- 
work, an individual is imagined to have latent FHT's 
Si,...,Sc for the C competing causes. The FHT S c 
is defined as in (1) with B c replacing B. The ob- 
served cause d and observed FHT Sd are then given 
by 

(3) Sd = min{S c , c = 1, . . . , C}. 

In other words, the observed FHT time and mode 
are those of the smallest latent FHT. 

As a practical matter, latent FHT's, other than 
the smallest, are generally unobservable, although 
there may be exceptions. For example, an engineer- 
ing system may consist of C components and fail 
whenever one of its components fails. If the sys- 
tem is repairable, then the failure times of all C 
components may eventually be observed because re- 
pair allows the system to continue functioning. It 
is questionable, however, if a repaired system re- 
ally continues on the same stochastic trajectory as 
before the repair. Nevertheless, latent FHT's have 
value as notional measures of survival increments 
that might be realized in an idealized world where 
selected causes of failure are removed. In terms of 
the definition in (3), an investigator might look at 
the impact on the observed FHT if some boundary 
subsets B c , associated with selected causes, are re- 
moved from the model formulation. This removal, in 
effect, eliminates the associated latent FHT's from 
(3). This exercise therefore simulates the elimina- 
tion of some causes of death (in medicine) or fail- 
ure (in engineering). An FHT model offers a clear 
conceptual structure for the status of an individual 
with respect to modes of failure other than the one 
observed. This structure is captured by the "dis- 
tance" that the threshold state X(Sd) lies from the 
boundary set B c for any cause c that differs from the 
observed cause d. The difference S c — Sd for c ^ d 
reflects the same distance on the survival time scale. 
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FHT models offer the possibility of inferences about 
this distance. 

To illustrate an FHT competing risks model by 
a concrete medical example, consider a multidimen- 
sional Wiener process of C dimensions with a bound- 
ary B c in each dimension. Each of the C dimensions 
defines a different cause of death c E {1, 2, . . . , C}. In 
such a model, one may make inferences about sec- 
ondary medical conditions that are not the primary 
cause of death. For example, in studies of occupa- 
tional exposure to diesel exhaust, workers may be 
found to have increased risks of death from lung 
disease, cardiovascular disease and other causes. It 
is desirable in such a context to have an FHT model 
that is capable of considering different causes of death 
simultaneously. A worker dying of lung cancer (the 
primary cause of death) may have advanced cardio- 
vascular disease, both of which are aggravated by 
exposure to an occupational hazard such as diesel 
exhaust. Then, an investigator's interest may lie in 
making inferences about the worker's cardiovascu- 
lar disease status at the time of death from lung 
cancer. We note that if the underlying multidimen- 
sional Wiener process is correlated, then the latent 
survival times for different causes of death will be 
dependent. 

5. CURE RATES 

We mentioned at the outset that some FHT mod- 
els may offer a positive probability of no FHT tak- 
ing place in finite time. Thus, for example, a med- 
ical treatment may offer a cure, some animals in a 
population may be immune to infection, some stock 
prices may never reach $100, and some marriages 
may never end in divorce. The fact that P(S = oo) > 
in some FHT models is closely related to compet- 
ing risks. Generally, if the FHT model takes account 
of all competing risks, then eventual failure from 
some cause is assured. If, however, the FHT model 
takes account of only one or a few competing risks, 
then there is a positive probability that the FHT will 
be infinite to accommodate those individuals who 
are not susceptible to the limited array of causes of 
failure that are considered in the model. To illustrate 
the natural way in which FHT models take account 
of a cure rate, consider a Wiener diffusion model 
with a fixed boundary at zero (the time axis). If the 
drift of the process is away from the boundary, that 
is, fj, > 0, then a finite FHT is not assured and, in 



particular, P(S < oo) = exp(— 2xQfi/a 2 ). Likewise, a 
gamma process with a cure rate might be defined as 

(4) X(t) = { X °' With P robabilit y 1 ~ P' 
' ^ ' \xq — Z(t), with probability p. 

Here the parameter p is a susceptibility fraction, 
with 0<p<l. As an example of this last model, 
a subject may have a malignant or benign form of a 
disease with probabilities p and 1 — p, respectively. 
The malignant form progresses monotonically to- 
ward a medical endpoint (e.g., death). 

6. COVARIATES AND LINK FUNCTIONS FOR 
THRESHOLD REGRESSION 

The parent process {X(t)} and boundary set B of 
the FHT model will both generally have parameters 
that depend on covariates that vary across individu- 
als. To illustrate, consider the Wiener process model 
in Example 3. The Wiener process has mean param- 
eter n and variance parameter a 2 and the boundary 
set has parameter xo, the initial process level. In 
threshold regression, these parameters will be con- 
nected to linear combinations of covariates using 
suitable regression link functions, as illustrated be- 
low for some parameter, say 6, 

(5) g e (e i )=z i (3. 

Here gg is the link function, the parameter 0{ is 
the value of the parameter 6 for individual i, Zj = 
(1, zn, . . . , Zik) is the covariate vector of individual 
i (with a leading unit to include an intercept term) 
and (3 is the associated vector of regression coef- 
ficients. The mathematical form of the link func- 
tion must be suited to the application. Generally, it 
will be chosen to map the parameter space into the 
real line. For example, a variance parameter such 
as <t may employ a logarithmic link function, that 
is, ln(c 2 ) = z(3. Likewise, the list of covariates and 
their mathematical forms in the regression function 
z/3 must be chosen appropriately, as is the case in a 
conventional regression analysis. 

Previous work that has considered regression struc- 
tures for FHT models includes Whitmore (1983), 
Whitmore, Crowder and Lawless (1998), Lee, 
DeGruttola and Schoenfeld (2000) and Lee et al. 
(2004). To illustrate one of these applications, Lee, 
DeGruttola and Schoenfeld (2000) use a bivariate 
Wiener diffusion process as the basis of a thresh- 
old regression model for the study of progression to 
death in AIDS, with CD4 cell count serving as a 
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marker process (marker processes are discussed in 
a later section). The initial health status and mean 
parameter of the parent process are made to de- 
pend on baseline covariates and treatment variables 
through log-linear and identity link functions, re- 
spectively. The mean and variance parameters of the 
marker process are also given a regression structure, 
with identity and log-linear link functions. Finally, 
the correlation parameter for the parent and marker 
processes uses a correlation transform as a link func- 
tion. 

Threshold regression raises some new issues for es- 
timation and inference in FHT models. Where FHT 
models are estimated only from censored survival 
data, parameter estimators may exhibit significant 
multicollinearity, especially with highly parameter- 
ized regression functions. This fact does not reflect 
any deficiency of the FHT model but, rather, reflects 
the limited information content of sample data in 
a rich modeling context. Reparameterization of the 
model can assist with computational problems that 
may arise from this multicollinearity but generally 
the condition is not sufficiently severe to prevent es- 
timates from being computed. The impact is primar- 
ily felt in the interpretability of the estimation re- 
sults. As with conventional regression, where regres- 
sion effects are highly collinear, it will be difficult to 
attribute the effect to a particular model compo- 
nent. For example, in threshold regression based on 
censored inverse Gaussian survival data, estimates 
of covariate effects of the initial value xq and mean 
parameter [i may be collinear because the mean sur- 
vival time depends on their ratio xq/\/j,\. Thus, the 
high correlation of their sampling errors can only 
be mitigated by having fine details for the disper- 
sion pattern of the survival data. Censoring or small 
sample sizes may mask those fine details and thus 
make estimation more difficult. 

7. RUNNING TIME VERSUS CALENDAR 
TIME 

In many applications of threshold regression, the 
natural time scale of the parent process is not calen- 
dar or clock time. For example, a mechanical compo- 
nent may wear according to the amount of its usage 
or liver disease may progress according to an individ- 
ual's cumulative consumption of alcohol. Mathemat- 
ical research on different time scales has been car- 
ried out by many researchers. Cox and Oakes (1984, 
Section 1.2, pages 3-4) pointed out that "often the 



'scale' for measuring time is clock time, although 
other possibilities certainly arise, such as the use of 
operating time of a system, mileage of a car, or some 
measure of cumulative load encountered." These ac- 
cumulation measures are increasing with calendar 
time and thus are alternative progression scales for 
the stochastic process. Such measures are given a 
variety of names, depending on the context, such 
as operational time, disease progression, or running 
time. We shall mainly use the last name here. If r(t) 
denotes the transformation of calendar time t to run- 
ning time r, with r(0) = 0, and {X(r)} is the par- 
ent process defined in terms of running time r, then 
the resulting process expressed in terms of calen- 
dar time is the subordinated process X*(t) = X[r(t)], 
where the asterisk identifies the subordinated pro- 
cess. Adaptations of FHT models to running time 
scales may be done in a variety of ways, as we de- 
scribe below. The variety includes both random and 
nonrandom transformations. We note that r(t) must 
be a monotonic transformation but its monotonic- 
ity need not be strictly increasing. Interesting ef- 
fects arise, for example, where r(t) is a function with 
jump discontinuities. 

1. Some applications require a monotonic mathe- 
matical transformation of the time scale. In these 
cases, r(t) is a deterministic function of calendar 
time t. A typical example from an engineering 
application is the strictly monotonic transforma- 
tion r = 1 — exp(— At 7 ) with A > and 7 > 0. See, 
for example, Carey and Koenig (1991), Whitmore 
(1995) and Whitmore and Schenkelberg (1997). 
The mathematical transformation may depend 
on covariates, as in Bagdonavicius and Nikulin 
(2001), where the running time scale forms part 
of an accelerated life model. 

2. Running time may also enter an FHT model us- 
ing a stochastic process for subordination. In this 
context, the parent process {X(r)} is directed by 
a second stochastic process {R(t)} having mono- 
tonic sample paths. In this context, we refer to 
{R(t)} as the directing process and the subordi- 
nated process takes the form {X*(t)} = {X[R(t}]}. 
Unlike a monotonic mathematical transformation, 
subordination with a stochastic process gives the 
transformation random properties that can greatly 
enrich the model. Lee and Whitmore (1993) ex- 
amine the connection between subordinated stochas- 
tic processes and running time. As a specific ex- 
ample of a subordinated process, one can con- 
sider a Poisson parent process that is directed by 
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a gamma process (which has monotonic sample 
paths). The result is a clustering Poisson process 
(Hougaard, Lee and Whitmore, 1997) in which 
an FHT can be triggered by the occurrence of a 
cluster of Poisson events. 
3. The running time scale may be a combination 
of different accumulation measures. For exam- 
ple, Oakes (1995) and Kordonsky and Gertsbakh 
(1997) look at multiple running time scales in sur- 
vival data analysis. Duchesne and Lawless (2000) 
and Duchesne and Rosenthal (2003) describe var- 
ious advances in running time models for survival 
data. The concept of collapsible time within the 
context of accelerated failure time models is cen- 
tral to this earlier work. The basic idea appears in 
various forms. For example, a composite running 
time might be defined by 

J 

(6) r{t)=Y j a j r ] {t), 

3=1 

where the rj(t) are different accumulation mea- 
sures that can advance degradation or disease 
progression and the ay are positive parameters 
that weight the contributions of the different mea- 
sures. One of the measures, say r\(t), may be cal- 
endar time itself so r\{t) = t. One ay parameter 
will need to be set to unity to give a well-defined 
scale. Typically, in this setup, composite running 
time has a fixed mathematical form for any given 
individual case but individuals will have different 
scales because the ry (t) vary randomly among in- 
dividuals. As a simple example of a composite 
running time, consider the mechanical aging of 
a motor vehicle which may be related to both 
the passage of calendar time r\{t) = t and accu- 
mulated mileage ^(i). In this case, (6) has two 
components, as follows: r{t) = t + aT2(i). Notice 
ot\ is set to 1 and «2 = ot. 

A practical case of the last kind of running time is 
illustrated in Lee et al. (2004) where railroad work- 
ers are employed in different types of jobs, indexed 
by j = 1, . . . , J, which have differential exposures 
to diesel exhaust, an occupational risk. The run- 
ning time (6) here is defined as a weighted sum of 
different exposure intervals. The quantity rj(t) is 
the time spent by the worker in job type j during 
time interval [0,t]. The ctj are positive weights that 
determine the rates at which the running time ad- 
vances per unit of calendar time spent in the differ- 
ent job types. One ay is set to unity as a numeraire, 



say ctj = 1. The ry (t) also must satisfy the account- 
ing constraint J2j=i r j(t) = Observe that (6) is a 
deterministic transformation for any given set of ex- 
posure intervals rj(t) but that these intervals vary 
randomly from one worker to another according to 
their individual work histories. 

Some processes may be defined in terms of run- 
ning time from the outset. For example, in a Bernoulli 
process or a Markov chain, the progress parameter 
represents the sequence of trials or steps of the pro- 
cesses, respectively. These parameters may already 
be seen as reflecting a kind of running time. The 
mapping of calendar time to this running time, as 
represented by r(t) = r, is already implicit in the 
discrete progress parameter of the process. 

The running time scale r(t) is included in the FHT 
model in order to make the model a more valid rep- 
resentation of reality. With a correct specification 
of running time, one would expect health status 
or component strength to decline steadily and pre- 
dictably against the scale that measures the accu- 
mulating "wear and tear" of running time. In other 
words, X(r) would retain very little or no inher- 
ent variability if r(t) could be chosen carefully. This 
situation describes an ideal that is unattainable in 
most practical applications of FHT models but is a 
target of model building. 

8. INCORPORATING MARKER PROCESSES 
IN THRESHOLD REGRESSION 

A marker process refers to an external process 
that covaries with the parent process. It assists in 
tracking progress of the parent process if the par- 
ent process is latent or only infrequently observed. 
In this way, the marker process forms a basis for 
predictive inference about the status of the parent 
process and its progress toward an FHT. Marker 
processes may also be of scientific interest in their 
own right. As markers of the parent process, they of- 
fer potential insights into the causal forces that are 
generating the movements of the parent process. Ex- 
amples of marker processes include CD4 cell count 
for AIDS, blood pressure for cardiovascular disease, 
personal medical cost for health status, input drive 
current for a laser, and ambient temperature for 
equipment. 

The basic analytical framework for a marker pro- 
cess conceives of a bivariate stochastic process 
{X(r),Y(r)} where the parent process {X(r)} is 
one component process and the marker process {y(r)} 
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is the other. Both are assumed to be one-dimensional 
for convenience of exposition. They are also both 
defined on the running time scale r. We discuss the 
implications of this last point shortly. Whitmore, 
Crowder and Lawless (1998) look at failure infer- 
ence based on a bivariate Wiener model in which 
failure is governed by the FHT of a latent degra- 
dation process while auxiliary readings are available 
from a correlated marker process. As noted earlier, 
Lee, DeGruttola and Schoenfeld (2000) apply this 
bivariate marker model to CD4 cell counts in the 
context of AIDS survival. 

An application may offer a variety of marker proces- 
ses, say, {Yfc(r), k = 1, . . . , K}, that may be of po- 
tential scientific value. They can be studied sepa- 
rately or combined into a composite marker process. 
For marker processes that involve measurements, 
the following additive form for the composite marker 
might be appropriate: 



K 



(7) 



Y{r) = lQ + Y,lkY k {i 



k=l 



The concept of a composite marker was first pro- 
posed by Whitmore, Crowder and Lawless (1998) 
in an engineering context. The aim in constructing 
the composite marker process is to find that lin- 
ear combination of the K candidate markers that 
has the largest predictive correlation or association 
with the parent process. The 7^ parameters define 
the linear combination and these generally must be 
estimated. The approach is reminiscent of regres- 
sion analysis with the composite marker serving as 
a regression function to predict the parent process 
{X(r)}. Here 70 serves as the intercept term of the 
regression relationship. If the composite marker can 
mimic the parent process perfectly, then {X(r)} and 
{y(r)} will be perfectly correlated. An exact model 
for the preceding setup is a (k + l)-variate Wiener 
diffusion process in which the parent process is one 
component and the k markers are the remaining 
components. The conditional process {X{r)\Y\{r) = 
yi(r), . . . ,Yk(t) = y_ft-(r)} then defines an exact lin- 
ear regression structure. Where one is dealing with 
a parent process or marker processes that are not 
measurement processes, such as Markov chains, the 
concept of a composite marker process must be re- 
defined in a suitable manner. 

The FHT modeling framework has evolved in the 
literature to encompass three major components as 
shown in Figure 1, namely, an FHT model (parent 
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the CT scale 





Fig. 1. This conceptual framework shows the connections 
between the parent process (often a latent process), running 
time (RT) and an external marker process that is correlated 
with the parent process. Time subordination links calendar 
time ( CT) to running time. The threshold regression struc- 
ture stands in the background and is not displayed explicitly 
in the figure. 

process and boundary set) that defines the relevant 
endpoint, a running time scale and a marker process. 
The threshold regression (TR) structure stands in 
the background of the schematic in Figure 1 and 
allows the parameters of the various components 
in the figure to depend on baseline and other co- 
variates. Although the schematic shows only one 
marker, it is clear that there may be many. The 
framework in Figure 1 has several noteworthy fea- 
tures. For example, if a marker process has mono- 
tonic sample paths, it may serve either 3jS 8, marker 
or as a running time r(t), as may be deemed appro- 
priate by the investigator. The framework reminds 
us that a marker process {Y(r)} should be defined 
on the running time scale r when its correlation with 
the parent process {X(r)} is being considered. Thus, 
for example, if r{t) measures an individual's cumu- 
lative exposure to a potential carcinogen at time t 
and the marker y is a serum measurement for the in- 
dividual on a cancer-specific antigen at time t, then 
the serum reading y should be recorded func- 
tion of cumulative exposure r. In other words, the 
progress parameter of the serum marker process is 
cumulative exposure r, not calendar time t. 

We have said that the parent process is generally 
latent. This feature definitely is common in medical 
applications where inherent health condition cannot 
be observed (and, indeed, may be deemed unmea- 
surable). Marker processes are surrogates for health 
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status, especially if they are highly correlated with 
the underlying medical condition. These markers may 
range from biomedical measurement processes, such 
as serum measurements, to more qualitative pro- 
cesses, such as periodic subjective evaluations of health 
status by a patient or caregiver. In engineering sys- 
tems, there will be contexts in which wear and tear, 
for example, can be observed and measured. In many 
physical settings, however, only surrogates for the 
system condition are available. For example, the drive 
current of a laser is a marker for its physical con- 
dition but not perfectly correlated with it. Similar 
comments can be made about social and economic 
systems. The parent processes that define the FHT 
(e.g., business failure) may only be imperfectly mon- 
itored by a marker process (e.g., accounting mea- 
sures of solvency). We also add that marker pro- 
cesses may be leading, lagging or coincident with 
respect to the parent process and their phase will 
be important in predictive inference for the parent 
process and its FHT. 

9. DATA STRUCTURES FOR THRESHOLD 
REGRESSION MODELS 

The data structures of threshold regression studies 
vary widely. To be specific, we look at the case in 
which longitudinal observations in calendar time are 
available for the parent process and the covariate 
vector. In this case, the data structure for a single 
individual can be summarized as follows: 

Time points: 

= t <t 1 <---<t m , 
Failure codes: 
(8) /o = 0, fi = 0, . . . , / m _i = 0, f m = or 1, 
Readings on parent process: 

Covariate vectors: 

z 0> z l > • • • ; z m- 

Each individual has observation vectors of the form 
(tj,fj,Xj,Zj), j = 0, l,...,m, where t = < h < 
• • • < t m . Here tj is the time of the jth observation, 
fj is an indicator variable for whether the time tj 
is an FHT, Xj is the state of the process at time tj 
and Zj is the covariate vector of the jth observation 
for the individual. 

The data structures may have a variety of special- 
ized features, as illustrated below. 

1. The data sets usually consist of a sample of indi- 
viduals, i = 1, . . . , n, with individual parent proces- 



ses {Xi(t)} and boundary sets . The individ- 
ual processes are often assumed to be mutually 
independent. 

2. Where there are competing modes of failure, then 
the cause of failure d will be recorded for each 
individual. 

3. The final observation time t m for an individual is 
a random stopping time if f m = 1. Thus, t m = S 
and x m = X(S) if f m = 1. Here X(S) £ B is the 
threshold state realized by the individual at the 
FHT. If f m = 0, then time t m is a right-censoring 
time for the FHT, that is, t m < S. If i m _i < S < 
t m , then the survival time S is interval censored. 

4. The data are longitudinal if there is more than 
one reading available for some individuals, that 
is, if to > 1 for some individuals. 

5. If the parent process is latent, then the data set 
will have no observations Xj , although there may 
still be readings on the covariate vectors Zj. 

6. If the data set consists only of a single time t and 
failure indicator / for each individual, then the 
data set constitutes censored survival data. With 
a baseline covariate vector zq available, the data 
provide a basis for censored survival threshold 
regression. 

7. Let X(tj) be abbreviated Xj for any individual. 
The reading Xj on the parent process, for j < to, 
is a realization of the conditional random vari- 
able Xj | S > tj . The conditioning event is that the 
process has reached state Xj at time tj without 
experiencing an FHT. 

8. Where {X(t)} is a Markov process (which is the 
most common type of model), we have for any 
individual that 

P(Xj = Xj\xj-i, ...,xq,S> tj) 

= P(Xj = Xj\Xj-i = Xj-\,S > tj) 

for j < m. 

In other words, the distribution of the next ob- 
servation Xj depends only on the value of the 
preceding observation Xj-\ and the fact that no 
FHT has yet occurred. The sample path by which 
Xj-i was attained is immaterial. 

Our discussion of data structures here has referred 
only to calendar time t without reference to run- 
ning time r. It also has not taken account of ob- 
servations that may be available on relevant marker 
processes. The discussion in Sections 7 and 8, how- 
ever, will make it clear what supplemental data are 
required when these model components are part of 
the threshold regression model. 
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10. PARAMETER ESTIMATION AND 
INFERENCE 

In applications to date, parameter estimation for 
FHT models and threshold regression have been heav- 
ily dominated by maximum likelihood methods. The 
reason is that the probabilistic specification of the 
parent stochastic process in FHT models is usually 
explicit and, hence, likelihood expressions follow as a 
matter of course. The optimization required by this 
estimation method may employ a variety of compu- 
tational techniques but gradient methods work very 
well. Extensions to Bayesian methods have been de- 
veloped in some cases. For example, Pettit and Young 
(1999) and Shubina (2005a, 2005b) have embedded 
the Wiener diffusion FHT model in a Bayesian frame- 
work. Lee, DeGruttola and Schoenfeld (2000) have 
developed some predictive inference results for this 
case, in conjunction with a marker process. Nothing 
seems to stand in the way of developing nonparamet- 
ric and semiparametric approaches for these models 
but these approaches have not yet been taken up in 
the literature. 

Case illustration. To illustrate the nature of in- 
ference for one of the simple threshold regression 
settings, we now set up the sample log-likelihood 
function for censored inverse Gaussian regression for 
a medical application like that found in Lee et al. 
(2004). We consider a latent health status process 
defined on a running time scale r. We let the par- 
ent process be a Wiener diffusion process. The FHT 
for such a process follows an inverse Gaussian distri- 
bution. The inverse Gaussian distribution depends 
on the mean and variance parameters of the un- 
derlying Wiener process (fi and a 2 ) and the initial 
health status level (xq). We let f{r\fi, a 2 , xq) and 
F(r\fi, a 2 , xq) denote the probability density func- 
tion (p.d.f.) and cumulative distribution function 
(c.d.f.) of the FHT distribution, both defined in terms 
of running time r. These functions have simple com- 
putational forms. For the case where the process be- 
gins at xq > and the boundary is the zero level, the 
p.d.f. for the first hitting time is given by 



P(X = oo) = 1 — exp(— 2xQfi/a 2 ). The c.d.f. corre- 
sponding to (9) is 

F(r\fi,a 2 ,x ) 

(fir + x ) 



f(r\fi,a 2 ,x ) 



X 



(9) 



V2-Ka 2 r 3 



exp 



(x + fir) 



2-i 



2a 2 i 



for — oo < fi < oo, a 2 > 0, xq > 0. 



If fi > 0, then the FHT is not certain to occur and 
the p.d.f. is improper. Specifically, in this case, 



(10) 



+ exp(-2x Q fi/a 2 )^ 



fir — Xq 



where <!>(•) is the c.d.f. of the standard normal dis- 
tribution. 

The health status process is latent here and, hence, 
can be given an arbitrary measurement unit. Thus, 
one parameter may be fixed. We set the variance pa- 
rameter a 2 to unity. Both parameters fi and xq are 
linked to k regression covariates that are represented 
by the row vector z = (1, z\, . . . , z k )- The leading 1 
in z allows for a constant term in the regression re- 
lationship. An identity function of the form 

fji = Z f3 = (3 + (3 1 z 1 H h f3 k z k 

is used to link the parameter fi to the covariates and 
a logarithmic function 

m(x ) = Z7 = 70 + 7i2i H h 7fcZfc 

is used to link the parameter xq to the covariates. 
Here f3 = (/3 , A, ■ • ■ , AO' and 7 = (70, 7i> • ■ • > 7fc)', 
where (3q and 70 are regression constants. Parame- 
ters of the running time scale, such asa = (a\ , . . . , a j) 
in (6), may also be linked to covariates using link 
functions of appropriate form. 

We now denote fi and xq for subject i by fi^> 
and Xq . We let r(0 denote the running time for sub- 
ject i. Time r(0 is the running time at the moment 
of death for a dying subject and a right-censored 
running time for the moment of death for a sur- 
viving subject. Hence, each dying subject i con- 
tributes probability density f(r^\fi^\xQ ) to the 
sample likelihood function, for i = 1, . . . ,n±, and each 
surviving subject i contributes survival probability 



W l,..(0 ^ 



F(rW\n 







1 - F(i 



(i)|.,(0 JO> 



I/' 







to the sample 
likelihood function, for i = m + 1, . . . , n\ + Uq. The 
sum n = n\ + riQ is the total number of subjects. 
The sample log-likelihood function to be maximized 
therefore has the form 

ni 



lnL( Q ,/3, 7 ) = E ln /0 



(0l„(0 JO 



/' 







'111 



ni+n 
i=n\+l 

Numerical gradient methods can be used to find 
maximum likelihood estimates for /3, 7 and a. 
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11. THRESHOLD REGRESSION FOR 
LONGITUDINAL DATA ANALYSIS 

Our discussion of data structures in Section 9 has 
anticipated that longitudinal data are gathered on 
the respective stochastic processes of individuals in 
some applications. Using our previous notation, we 
now let {Aj} denote the longitudinal observation 

process, defined on the time points tj,j = 0, 1, If 

the individual survives beyond time tj, then the fail- 
ure code fj = and Aj = {S > tj,Xj,Zj} for j < m. 
If the individual fails in the final interval (t m -i,t m ], 
then f m = 1 and A m = {S G (t m -i,t m ],x m G B}. As 
defined earlier, S is the stopping time for the longi- 
tudinal observation process. We note that z m is not 
defined when the individual has failed and, hence, is 
dropped from the expression for A m . Moreover, the 
final reading x m for the parent process lies inside 
the boundary set B when the individual has failed. 

Longitudinal data of this kind pose an interesting 
challenge for first-hitting-time models, as for most 
time-to-event models. Lu (1995) considers the prob- 
lem for the basic Wiener model where longitudinal 
observations are made on the process {X(t)} up to 
the hitting or censoring time, as the case may be. 
She formulates the likelihood function and computes 
exact maximum likelihood estimates. The method- 
ology is somewhat intricate but manageable. Lee, 
DeGruttola and Schoenfeld (2000) consider the is- 
sue of modeling longitudinal data for a bivariate 
Wiener model representing a latent health status 
process and a correlated marker process. These au- 
thors mention an interesting approach to handling 
longitudinal data which they anticipated would be 
technically satisfactory and practical to implement. 
Their suggested approach, however, is not elabo- 
rated in their article, so we sketch one direction of 
development below but leave a full exploration of 
the approach as an open research question. We refer 
to this method as an uncoupling procedure because 
it effectively unlinks the longitudinal observations 
into a set of independent conditional observations. 

With the preceding notation, the probability of 
observing the longitudinal data record of an indi- 
vidual can be expanded as a product of conditional 
probabilities as 

P(A m ,A m - 1 ,...,A 1 ,A ) 

(12) 

m 

= P(A )l[P(A j \A j - 1 ,...,A ). 



Now we come to the crucial assumption. If it can be 
assumed that {Aj,j = 0, 1, . . .} is a Markov process 
with initial state Aq, then (12) can be simplified as 

P(A m ,A m _ 1 ,...,A 1 ,A ) 

(13) 

m 

= P(A Q )l[P(A j \A j - 1 ). 

i=i 

In other words, the probability of observing Aj de- 
pends only on its preceding state and not on 
the earlier history of the observation process. The 
explicit forms of the probability elements on the 
right-hand side of (13) are 

P{Aj\Aj-{) 

(14) = P(S > tj,Xj,Zj\S > tj_i,Xj_i,Zj_i) 

if fj =0,j <m, 

P\A m \A m — i) 

(15) =P{Se{t m - U t m ],x m eB\ 

S > t m — 1; X m — i, Z m ~ i) if f m = 1. 

If no observations are available on the parent pro- 
cess, then Xj is dropped from the Aj notation, giving 
Aj = {S > tj,Zj} if fj = 0, j < to, and A m = {S G 
(i m _i,t m ]} if f m = 1. Again, invoking the Markov 
assumption for the observation process, (14) and 

(15) take the revised forms 

P(Aj\Aj- X ) 

(16) =P(S>tj,Zj\S>tj- 1> z j - 1 ) 

if fj = for j < m, 

P(A m \A m — i) 

(17) =P(Se(Vi 

if f m = 1. 

Statement (13) is the theoretical justification for the 
uncoupling procedure. Neither this theoretical devel- 
opment nor issues of practical implementation of the 
procedure were taken up by Lee, DeGruttola and 
Schoenfeld (2000). As already noted, the procedure 
remains an open topic for future research. 

12. MODEL VALIDATION, DIAGNOSTICS 
AND REMEDIES 

Although procedures for model validation, diag- 
nostics and remedies are not as well developed for 
threshold regression as for conventional regression 
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models for survival data, a number of techniques 
have been proposed and applied successfully in ear- 
lier FHT investigations. For example, procedures are 
available for checking the assumptions of the TR re- 
gression model having a Wiener process and inverse 
Gaussian FHT, both with and without associated 
marker processes. Lee, DeGruttola and Schoenfeld 
(2000) present some procedures for this TR model 
and demonstrate the techniques using a medical case 
application. Lee and Whitmore (2002) present a 
larger suite of techniques for checking assumptions 
of this model and also discuss a number of reme- 
dies that might be used where assumptions do not 
hold. Lee et al. (2004) also discuss validation for an 
extension of this same model in which the calen- 
dar time scale is replaced by a job-exposure disease 
progression scale. One of the proposed validation 
procedures relies on the fact that the inverse Gaus- 
sian (IG) distribution is the first-stopping-time dis- 
tribution of a Wiener process. Hence, comparisons 
can be made between Kaplan-Meier (KM) survival 
curves and the IG survival curves implied by the 
model (for different covariate subgroups). Applica- 
tions with longitudinal observations on the parent 
process or marker measurements offer even more 
data for model validation. The previous work also 
points out the importance of having subject-matter 
specialists understand the model features and com- 
pare them with the fundamental physical processes 
at play. For example, the concept of an FHT is one 
feature whose mechanism is found frequently in na- 
ture, is easily understood by scientists and can be 
checked against their scientific understanding of the 
application context. 

13. SOME OPEN RESEARCH PROBLEMS 

Many interesting aspects of threshold regression 
require further study. We noted earlier that multi- 
collinearity of parameter estimates can be a prac- 
tical issue. It remains to be seen which parame- 
terizations of threshold regression models tend to 
have relatively independent estimation errors. Mul- 
ticollinearity within regression functions will tend to 
show itself in familiar ways and will likely be dealt 
with by conventional remedies. 

The Cox proportional hazards regression model 
is widely used for survival data analysis. Thresh- 
old regression models do not generally possess the 
proportional hazards feature for different configu- 
rations of covariates. A useful research contribution 



would be made by comparing and contrasting the re- 
sults of Cox regression and TR in the same context. 
Some public sets of survival data that are scientifi- 
cally important and have a plausible FHT interpre- 
tation might very well be reanalyzed to see if the key 
research conclusions are materially affected when a 
TR model is used in place of a more conventional 
technique. 

Both parent and marker processes may be sub- 
ject to measurement error. For example, blood pres- 
sure is known to be measured with error. Whitmore 
(1995), for example, studied a Wiener diffusion FHT 
model with measurement error. The true state of 
a process is also often randomly masked. The in- 
corporation of measurement or masking errors in 
TR models, where these extensions are motivated 
by significant applications, would represent a useful 
research extension. 

The identification of individual marker processes 
and the construction of composite marker processes 
to track or mimic a latent parent process are chal- 
lenging subjects that need further theoretical work 
and more experience with real applications. The chal- 
lenge will be especially great where the marker pro- 
cesses and latent processes are from different classes 
of processes. A related open research issue concerns 
the investigation of whether markers are leading, 
lagging or coincident with the parent process. 

Nonparametric, semiparametric and other robust 
estimation methods seem to have much to contribute 
to the successful application of threshold regression. 
Quasi- likelihood methods and generalized estimat- 
ing equations may offer feasible approaches. As thresh- 
old regression estimation in a general setting in- 
volves parameter estimation for the boundary set, 
the parent process and the running time scale, it 
is conceivable that a blend of nonparametric and 
parametric methods may be effective in some appli- 
cations. For example, nonparametric estimation of 
running time parameters might be combined with 
parametric estimation of the parent process. 

Our discussion of the analysis of longitudinal data 
in the context of threshold regression has already 
pointed out that a full theoretical development and 
justification of the uncoupling method remains an 
open research issue. In the same vein, practical ex- 
perience with this method or other methods for han- 
dling longitudinal data in threshold regression will 
be valuable contributions. 

Much remains to be done on model validation and 
diagnostic techniques in the context of threshold re- 
gression. These tools are likely to be developed as 
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threshold regression is applied in a broader range of 
practical cases. The earlier work on model valida- 
tion has been largely restricted to the Wiener FHT 
model and thus extensions to other FHT models 
need attention. For example, comparisons of Kaplan- 
Meyer (KM) survival curves with fitted TR survival 
curves, both defined on running time scales, will re- 
quire new methods that take account of the fact that 
the running time scale is itself fitted by a statisti- 
cal model. As another example, TR models assume 
that particular functions link the model parameters 
to the regression covariates. Both the forms of the 
link functions and the adequacy of the regression 
functions must be validated. Whether the correct di- 
recting process has been chosen is also a feature that 
must be checked by model validation techniques. Al- 
though model validity is likely to be established by 
standard techniques (such as cross-validation), new 
techniques and modifications of conventional meth- 
ods will surely be needed. In addition to using sta- 
tistical methods for model verification, it is desirable 
to work closely with subject-matter specialists to en- 
sure that the FHT models have realistic features and 
that the findings emerging from the analysis make 
practical sense. 

The last sentence of the preceding paragraph hints 
at the largest open research question. Threshold re- 
gression will prove itself through beneficial practical 
application. With exploration of fresh application 
areas will come ideas for better methods and mod- 
els for this new type of regression approach. 
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